###############################
# Panel Match: Crime Concerns
###############################

rm(list=ls(all=TRUE)) 

library(PanelMatch)
set.seed(1111)

#sink("15_panel_match.txt")

# load data
load("final_county_2023april6.RData")
names(d_county)

# generate dataset
d = data.frame(d_county$survey,
               d_county$county,
               d_county$t_stagg_haiti,
               d_county$t_stagg_all,
               d_county$crime_concerns_s,
               d_county$crime_rates_s,
               d_county$turnout_2013,
               d_county$income_2003,
               d_county$population_2011,
               d_county$right_voteshare_2013)
colnames(d) = c("survey",
                "county",
                "t_stagg_haiti",
                "t_stagg_all",
                "crime_concerns_s",
                "crime_rates_s",
                "turnout_2013",
                "income_2003",
                "population_2011",
                "right_voteshare_2013")
names(d)

# check treatment
table(d$t_stagg_haiti)

#############################################
# Matching Mahalanobis
#############################################

# Concerns Haiti
pm.results1 <- PanelMatch(time.id = "survey", 
                          unit.id = "county", 
                          treatment = "t_stagg_haiti", 
                          outcome.var = "crime_concerns_s",
                          refinement.method = "mahalanobis", 
                          covs.formula = ~turnout_2013 + income_2003 + population_2011 + right_voteshare_2013, 
                          match.missing = T, 
                          lag = 7, 
                          size.match = 1, 
                          qoi = "att" ,
                          forbid.treatment.reversal = FALSE,
                          data = d)

pe.results1 <- PanelEstimate(sets = pm.results1, data = d)
results1 = summary(pe.results1)[1]

# Concerns all
pm.results2 <- PanelMatch(time.id = "survey", 
                          unit.id = "county", 
                          treatment = "t_stagg_all", 
                          outcome.var = "crime_concerns_s",
                          refinement.method = "mahalanobis", 
                          covs.formula = ~turnout_2013 + income_2003 + population_2011 + right_voteshare_2013, 
                          match.missing = T, 
                          lag = 7, 
                          size.match = 1, 
                          qoi = "att" ,
                          forbid.treatment.reversal = FALSE,
                          data = d)

pe.results2 <- PanelEstimate(sets = pm.results2, data = d)
results2 = summary(pe.results2)[1]

# Rates Haiti
pm.results3 <- PanelMatch(time.id = "survey", 
                          unit.id = "county", 
                          treatment = "t_stagg_haiti", 
                          outcome.var = "crime_rates_s",
                          refinement.method = "mahalanobis", 
                          covs.formula = ~turnout_2013 + income_2003 + population_2011 + right_voteshare_2013, 
                          match.missing = T, 
                          lag = 7, 
                          size.match = 1, 
                          qoi = "att" ,
                          forbid.treatment.reversal = FALSE,
                          data = d)

pe.results3 <- PanelEstimate(sets = pm.results3, data = d)
results3 = summary(pe.results3)[1]

# Rates all
pm.results4 <- PanelMatch(time.id = "survey", 
                          unit.id = "county", 
                          treatment = "t_stagg_all", 
                          outcome.var = "crime_rates_s",
                          refinement.method = "mahalanobis", 
                          covs.formula = ~turnout_2013 + income_2003 + population_2011 + right_voteshare_2013, 
                          match.missing = T, 
                          lag = 7, 
                          size.match = 1, 
                          qoi = "att" ,
                          forbid.treatment.reversal = FALSE,
                          data = d)

pe.results4 <- PanelEstimate(sets = pm.results4, data = d)
results4 = summary(pe.results4)[1]

###############################################
# Table
###############################################

d1 = data.frame(results1)
d2 = data.frame(results2)
d3 = data.frame(results3)
d4 = data.frame(results4)

pe1 = round(d1$summary.estimate,3)
c1_lower = round(d1$summary.2.5.,3)
c1_upper = round(d1$summary.97.5.,3)

pe2 = round(d2$summary.estimate,3)
c2_lower = round(d2$summary.2.5.,3)
c2_upper = round(d2$summary.97.5.,3)

pe3 = round(d3$summary.estimate,3)
c3_lower = round(d3$summary.2.5.,3)
c3_upper = round(d3$summary.97.5.,3)

pe4 = round(d4$summary.estimate,3)
c4_lower = round(d4$summary.2.5.,3)
c4_upper = round(d4$summary.97.5.,3)

# Table a7 
columns = c("visas from Haiti","visas from all countries")
point_estimate7 = c(pe1,pe2)
conf_inter_lower7 = c(c1_lower,c2_lower)
conf_inter_upper7 = c(c1_upper,c2_upper)
table7 = cbind(columns,point_estimate7,conf_inter_lower7,conf_inter_upper7)
table7

write.table(table7, "tableA7.txt", sep="\t")

# Table a8
point_estimate8 = c(pe3,pe4)
conf_inter_lower8 = c(c3_lower,c4_lower)
conf_inter_upper8 = c(c3_upper,c4_upper)
table8 = cbind(columns,point_estimate8,conf_inter_lower8,conf_inter_upper8)
table8

write.table(table8, "tableA8.txt", sep="\t")

sink()

